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1. Introduction 

Turbulence is often described as the last great unsolved problem of classical physics [1-3]. 
However, it is not easy to state what would constitute a solution of the turbulence prob- 
lem. This is principally because turbulence is not one problem but a collection of several 
important problems: These include the characterisation and control of turbulent flows, 
both subsonic and supersonic, of interest to engineers such as flows in pipes or over cars 
and aeroplanes [4,5]. Mathematical questions in this area are concerned with develop- 
ing proofs of the smoothness, or lack thereof, of solutions of the Navier-Stokes and re- 
lated equations [6-10]. Turbulence also provides a variety of challenges for fluid dy- 
namicists [5,1 1-13], astrophysicists [14-17], geophysicists [18,19], climate scientists [20], 
plasma physicists [15-17,21,22], and statistical physicists [23-32]. In this brief overview, 
written primarily for physicists who are not experts in turbulence, we concentrate on some 
recent advances in the statistical characterisation of fluid turbulence [33] in three dimen- 
sions, the turbulence of passive scalars such as pollutants [34], two-dimensional turbu- 
lence in thin films or soap films [35,36], turbulence in the Burgers equation [37-39], and 
fluid turbulence with polymer additives [40^-2]; in most of this paper we restrict our- 
selves to homogeneous, isotropic turbulence [33,43,44]; and we highlight some similar- 
ities between the statistical properties of systems at a critical point and those of turbu- 
lent fluids [31,45,46]. Several important problems that we do not attempt to cover include 
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Rayleigh-Benard turbulence [47], superfluid turbulence [3,48], magnetohydrodyanmic tur- 
bulence [15,17,21,22], the behaviour of inertial particles in turbulent flows [49], the transi- 
tion to turbulence in different experimental situations [50,51], and boundary-layer [52,53] 
and wall-bounded [54] turbulence. 

This paper is organised as follows: Section 2 gives an overview of some of the experi- 
ments of relevance to our discussion here. In Section 3 we introduce the equations that we 
consider. Section 4 is devoted to a summary of phenomenological approaches that have 
been developed, since the pioneering studies of Richardson [55] and Kolmogorov [56], in 
1941 (K41), to understand the behaviour of velocity and other structure functions in in- 
ertial ranges. Section 5 introduces the ideas of multiscaling that have been developed to 
understand deviations from the predictions of K41-type phenomenology. Section 6 con- 
tains illustrative direct numerical simulations; it consists of five subsections devoted to 
(a) three-dimensional fluid turbulence, (b) shell models, (c) two-dimensional turbulence in 
soap films, (d) turbulence in the one-dimensional Burgers equation, and (e) fluid turbulence 
with polymer additives. Section 7 contains concluding remarks. 

2. Experimental Overview 

Turbulent flows abound in nature. They include the flow of water in a garden pipe or in 
rapids, the flow of air over moving cars or aeroplanes, jets that are formed when a fluid is 
forced through an orifice, the turbulent advection of pollutants such as ash from a volcanic 
eruption, terrestrial and Jovian storms, turbulent convection in the sun, and turbulent shear 
flows in the arms of spiral galaxies. A wide variety of experimental studies have been 
carried out to understand the properties of such turbulent flows; we concentrate on those 
that are designed to elucidate the statistical properties of turbulence, especially turbulence 
that is, at small spatial scales and far away from boundaries, homogeneous and isotropic. 
Most of our discussion will be devoted to incompressible flows, i.e., low-Mach-number 
cases in which the fluid velocity is much less than the velocity of sound in the fluid. 

In laboratories such turbulence is generated in many different ways. A common method 
uses a grid in a wind tunnel [57]; the flow downstream from this grid is homogeneous and 
isotropic, to a good approximation. Another technique use the von Karman swirling flow, 
i.e., flow generated in a fluid contained in a cylindrical tank with two coaxial, counterrotat- 
ing discs at its ends [58-60] ; in the middle of the tank, far away from the discs, the turbulent 
flow is approximately homogeneous and isotropic. Electromagnetically forced thin films 
and soap films [1,35,36] have yielded very useful results for two-dimensional turbulence. 
Turbulence data can also be obtained from atmospheric boundary layers [61-64], oceanic 
flows [65], and astrophysical measurements [14]; experimental conditions cannot be con- 
trolled as carefully in such natural settings as they can be in a laboratory, but a far greater 
range of length scales can be probed than is possible in laboratory experiments. 

Traditionally, experiments have measured the velocity u(x, t) at a single point x at var- 
ious times t by using hot-wire anemometers; these anemometers can have limitations in 
(a) the number of components of the velocity that can be measured and (b) the spatial and 
temporal resolutions that can be obtained [66,67]. Such measurements yield a time series 
for the velocity; if the mean flow velocity U >> u rms , the root-mean-square fluctuations 
of the velocity, then Taylor's frozen-flow hypothesis [5,33] can be used to relate temporal 
separations St to spatial separations Sr, along the mean flow direction via Sr = USt. The 
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Reynolds number Re = UL/v, where U and L are typical velocity and length scales in 
the flow and v is the kinematic viscosity, is a convenient dimensionless control parameter; 
at low Re flows are laminar; as it increases increases there is a transition to turbulence 
often via a variety of instabilities [50] that we will not cover here; and at large Re fully 
developed turbulence sets in. To compare different flows it is often useful to employ the 
Taylor-microscale Reynolds number Re\ = u rms \/v, where the Taylor microscale A can 
be obtained from the energy spectrum as described below (Sec. 6.3). 

Refinements in hot-wire anemometry [63,68] and flow visualisation techniques such as 
laser-doppler velocimetry (LDV) [66], particle-image velocimetry (PIV) [66,67], particle- 
tracking velocimetry (PTV) [66,67], tomographic PIV [69], holographic PIV [70], and 
digital holographic microscopy [71] have made it possible to obtain reliable measurements 
of the Eulerian velocity u(x, t) (see Sec. 3) in a turbulent flow. In the simplest forms of 
anemometry a time series of the velocity is obtained at a given point in space; in PIV two 
components of the velocity field can be obtained in a sheet at a given time; holographic 
PIV can yield all components of the velocity field in a volume. Components of the velocity 
derivative tensor = djUi can also be obtained [63] and thence quantities such as the 
energy dissipation rate per unit mass per unit volume e = —v^Yln j(diUj + djUi) 2 , the 
vorticity uu = V x u, and components of the rate of strain tensor = (diUj + djUi)/2, 
where the subscripts i and j are Cartesian indices. A discussion of the subtleties and lim- 
itations of these measurement techniques lies beyond the scope of our overview; we refer 
the reader to Refs. [63,66,67] for details. Significant progress has also been made over the 
past decade in the measurement of Lagrangian trajectories (see Sec. 3) of tracer particles 
in turbulent flows [58,59]. Given such measurements, experimentalists can obtain several 
properties of turbulent flows. We give illustrative examples of the types of properties we 
consider. 

Flow-visualisation methods often display large-scale coherent structures in turbulent 
flows. Examples of such structures plumes in Rayleigh-Bernard convection [72], struc- 
tures behind a splitter plate [73], and large vortical structures in two-dimensional or strat- 
ified flows [1,35,36]. In three-dimensional flows, as we will see in greater detail below, 
energy that is pumped into the flow at the injection scale L cascades, as first suggested 
by Richardson [55], from large-scale eddies to small-scale ones till it is eventually dis- 
sipated around and beyond the dissipation scale By contrast, two-dimensional turbu- 
lence [35,36,74,75] displays a dual cascade: there is an inverse cascade of energy from the 
scale at which it is pumped into the system to large length scales and a direct cascade of 
enstrophy ft = (^uj 2 ) to small length scales. The inverse cascade of energy is associated 
with the formation of a few large vortices; in practical realisations the sizes of such vor- 
tices are controlled finally by Ekman friction that is induced, e.g., by air drag in soap-film 
turbulence. 

Measurements of the vorticity u> in highly turbulent flows show that regions of large w 
are organised into slender tubes. The first experimental evidence for this was obtained by 
seeding the flow with bubbles that moved preferentially to regions of low pressure [76] that 
are associated with large-w regimes. For recent experiments on vortex tubes we refer the 
reader to Ref. [77]. 

The time series of the fluid velocity at a given point x shows strong fluctuations. It 
is natural, therefore, to inquire into the statistical properties of turbulent flows. From the 
Eulerian velocity u(x, t) and its derivatives we can obtain one-point statistics, such as 
probability distribution functions (PDFs) of the velocity and its derivatives. Velocity PDFs 
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are found to be close to Gaussian distributions. However, PDFs of lo 2 and velocity deriva- 
tives show significant non-Gaussian tails; for a recent study, which contains references to 
earlier work, see Ref. [63]. The PDF of e is non-Gaussian too and the time series of e 
is highly intermittent [78]; furthermore, in the limit Re — ► oo, i.e., v — > 0, the energy 
dissipation rate per unit volume e approaches a positive constant value (see, e.g., Fig. 2 of 
Ref. [79]), a result referred to as a dissipative anomaly or the zeroth law of turbulence. 

Various statistical properties of the rate-of-strain tensor, with components s%j, have been 
measured [63]. The eigenvalues Ai, A 2 , and A3, with Ai > A 2 > A 3 , of this tensor must 
satisfy Ai + A 2 + A3 = 0, with Ai > and A 2 < 0, in an incompressible flow. The sign of 
A 2 cannot be determined by this condition but its PDF shows that, in turbulent flows, A 2 has 
a small, positive mean value [80]; and the PDFs of cos(o> • e,), where is the normalised 
eigenvector corresponding to Ai, show that there is a preferential alignment [63] of 10 and 
e 2 . Joint PDFs can be measured too with good accuracy. An example of recent interest is a 
tear-drop feature observed in contour plots of the joint PDF of, respectively, the second and 
third invariants, Q = —tr{A 2 )/2 and R = — tr(A 3 )/3 of the velocity gradient tensor Aij 
(see Fig. 1 1 of Ref. [63]); we display such a plot in Sec. 6 that deals with direct numerical 
simulations. 

Two-point statistics are characterised conventionally by studying the equal-time, order- 
p, longitudinal velocity structure function 

S p (r) = ([(«(x + r)-u(x))-(r/rH, (1) 

where the angular brackets indicate a time average over the nonequilibrium statisti- 
cal steady state that we obtain in forced turbulence (decaying turbulence is discussed 
in Sec. 6.2). Experiments [33,81] show that, for separations r in the inertial range 
rj d « r « L, 

S p (r)~r& (2) 

with exponents ( p that deviate significantly from the simple scaling prediction [56] C, p A1 = 
p/3, especially for p > 3, where ( p < Cp 41 - This prediction, made by Kolmogorov in 
1941 (hence the abbreviation K41), is discussed in Sec. 4 below; the deviations from this 
simple scaling prediction are referred to as multiscaling (Sec. 5) and they are associated 
with the intermittency of e mentioned above. We mention, in passing, that the log-Poisson 
model due to She and Leveque provides a good parametrisation of the plot of ( p versus p 
[82]. 

The second-order structure function 5 2 (r) can be related easily by Fourier transforma- 
tion to the the energy spectrum E(k) = 4nk 2 (|u(fc) | 2 ), where the tilde denotes the Fourier 
transform, k = |k|, k is the wave vector, we assume that the turbulence is homogeneous 
and isotropic, and, for specificity, we give the formula for the three-dimensional case. 
Since Cf 41 = 2/3, the K41 prediction is 

E Kil (k) ~ fc- 5 / 3 , (3) 

a result that is in good agreement with a wide range of experiments [see, e.g., 
Refs. [33,83]]. 

The structure functions S p (r) are the moments of the PDFs of the longitudinal velocity 
increments 5u\\ = [(w(x + r) — u(x)) • (r/r)]. [In the argument of S p we use r instead of 
r when we consider homogeneous, isotropic turbulence.] These PDFs have been measured 
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directly [84] and they show non-Gaussian tails; as r decreases, the deviations of these 
PDFs from Gaussian distributions increases. 

We now present a few examples of recent Lagrangian measurements [58,59] that have 
been designed to track tracer particles in, e.g., the von Karman flow at large Reynolds 
numbers. By employing state-of-the-art measurement techniques, such as silicon strip de- 
tectors [59], used in high-energy-physics experiments, or acoustic-doppler methods [58], 
these experiments have been able to attain high spatial resolution and high sampling rates 
and have, therefore, been enable to obtain good data for acceleration statistics of La- 
grangian particles and the analogues of velocity structure functions for them. 

These experiments [59] find, for 500 < Re\ < 970, consistency with the Heisenberg- 
Yaglom scaling form of the acceleration variance, i.e., 

(a iaj ) ~ eWW-Wfiij, (4) 

where ai is the Cartesian component i of the acceleration. Furthermore, there are indica- 
tions of strong intermittency effects in the acceleration of particles and anisotropy effects 
are present even at very large Re\. 

Order-p Lagrangian velocity structure functions are defined along a Lagrangian trajec- 
tory as 

St P (.r) = ([vHt + r)-vf'(t)n (5) 

where the superscript L denotes Lagrangian and the subscript i the Cartesian component. 
If the time lag r lies in the temporal analogue of the inertial range, i.e., t v <C t <C T l , 
where t v is the viscous dissipation time scale and Tl is the time associated with the scale 
L at which energy is injected into the system, then it is expected that 

S&,(T)~rtf-. (6) 

The analogue of the dimensional K41 prediction is (^p K41 — p/2; experiments and simu- 
lations [60] indicate that there are corrections to this simple dimensional prediction. 

The best laboratory realisations of two-dimensional turbulence are (a) a thin layer of a 
conducting fluid excited by magnetic fields, varying both in space and time and applied 
perpendicular to the layer [85], and (b) soap films [86] in which turbulence can be gener- 
ated either by electromagnetic forcing or by the introduction of a comb, which plays the 
role of a grid, in a rapidly flowing soap film. In the range of parameters used in typical 
experimental studies [1,35,36,87] both these systems can be described quite well [88,89] 
by the 2D Navier Stokes equation (see Sec. 3) with an additional Ekman-friction term, 
induced typically by air drag; however, in some cases we must also account for corrections 
arising from fluctuations of the film thickness, compressibility effects, and the Marangoni 
effect. Measurement techniques are similar to those employed to study three-dimensional 
turbulence [1,35,36]. Two-dimensional analogues of the PDFs described above for 3D tur- 
bulence have been measured [see, e.g., Refs. [87]]; we will touch on these briefly when 
we discuss numerical simulations of 2D turbulence in Sec. 6.3. Velocity and vorticity 
structure functions can be measured as in 3D turbulence; however, inertial ranges associ- 
ated with inverse and forward cascades must be distinguished; the former shows simple 
scaling with an energy spectrum E(k) ~ /c~ 5 / 3 whereas the latter has an energy spectrum 
E(k) ~ k-^+ s \ with S = if there is no Ekman friction and 6 > otherwise. In the 
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forward cascade velocity structure functions show simple scaling [87]; we are not aware 
of experimental measurements of vorticity structure functions (we will discuss these in the 
context of numerical simulations in Sec. 6.3). 

We end this Section with a brief discussion of one example of turbulence in a non- 
Newtonian setting, namely, fluid flow in the presence of polymer additives. There are 
two dimensionless control parameters in this case: Re and the Weissenberg number We, 
which is a ratio of the polymer-relaxation time and a typical shearing time in the flow 
(some studies [41] use a similar dimensionless parameter called the Deborah number De). 
Dramatically different behaviours arise depending on the values of these parameters. 

In the absence of polymers the flow is laminar at low Re; however, the addition of 
small amounts of high-molecular- weight polymers can induce elastic turbulence [90], i.e., 
a mixing flow that is like turbulence and in which the drag increases with increasing We. 
We will not discuss elastic turbulence in detail here; we refer the reader to Ref. [90] for an 
overview of experiments and to Ref. [91] for representative numerical simulations. 

If, instead, the flow is turbulent in the absence of polymers, i.e., we consider large-i?e 
flows, then the addition of polymers leads to the dramatic phenomenon of drag reduction 
that has been known since 1949 [92]; it has obvious and important industrial applications 
[40,41,93-95]. Normally drag reduction is discussed in the context of pipe or channel 
flows: on the addition of polymers to turbulent flow in a pipe, the pressure difference 
required to maintain a given volumetric flow rate decreases, i.e., the drag is reduced and 
a percentage drag reduction can be obtained from the percentage reduction in the pressure 
difference. For a recent discussion of drag reduction in pipe or channel flows we refer the 
reader to Ref. [41]. Here we concentrate on other phenomena that are associated with the 
addition of polymers to turbulent flows that are homogeneous and isotropic. In particular, 
experiments [93] show that the polymers lead to a suppression of small-scale structures 
and important modifications in the second-order structure function [96]. We will return to 
an examination of such phenomena when we discuss direct numerical simulations in Sec. 
6.5. 



3. Models 

Before we discuss advances in the statistical characterization of turbulence, we provide a 
brief description to the models we consider. We start with the basic equations of hydrody- 
namics, in three and two dimensions, that are central to studies of turbulence. We also give 
introductory overviews of the Burgers equation in one dimension, the advection-diffusion 
equation for passive scalars, and the coupled NS and finitely extensible nonlinear elastic 
Peterlin (FENE-P) equations for polymers in a fluid. We end this Section with a description 
of shell models that are often used as highly simplified models for homogeneous, isotropic 
turbulence. 

At low Mach numbers, fluid flows are governed by the Navier-Stokes (NS) Eq. (7) 
augmented by the incompressibility condition 

d t u + (u.V)u = -Vp + vV 2 u + f , 

V • u = 0, (7) 

where we use units in which the density p = 1 , the Eulerian velocity at point r and time t is 
u(r, t), the external body force per unit volume is f , and v is the kinematic viscosity. The 
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pressure p can be eliminated by using the incompressibility condition [5,33,43] and it can 
then be obtained from the Poisson equation \7 2 p = —dij(uiUj). In the unforced, inviscid 
case, the momentum, the kinetic energy, and the helicity H = J drui • u/2 are conserved; 
here w = V x u is the vorticity. The Reynolds number Re = LV/v, where L and V are 
characteristic length and velocity scales, is a convenient dimensionless control parameter: 
The flow is laminar at low Re and irregular, and eventually turbulent, as Re is increased. 
In the vorticity formulation the NS equation 7 becomes 

d t uj = V x u x u + vV 2 uj + V x f ; (8) 

the pressure is eliminated naturally here. This formulation is particularly useful is two 
dimensions since lo is a pseudo-scalar in this case. Specifically, in two dimensions, the NS 
equation can be written in terms of w and the stream function tp: 

d t uj - J(ip, uj) = vS7 2 lu + oieu + /; 

Jty,U) = {d X ^)(dyUj) - {d X Lj)(dylP). (9) 

Here cxe is the coefficient of the air-drag-induced Ekman-friction term. The incompress- 
ibility constraint 

8 X U X + 8yUy = (10) 

ensures that the velocity is uniquely determined by ip via 

u = (-d y i>,d x ip). (11) 

In the inviscid, unforced case we have more conserved quantities in two dimensions than 
in three; the additional conserved quantities are for all powers n, the first of which 

is the mean enstrophy, il = {\u 2 ). 

In one dimension (ID) the incompressibility constraint leads to trivial velocity fields. 
It is fruitful, however, to consider the Burgers equation [37], which is the NS equation 
without pressure and the incompressibility constraint. This has been studied in great detail 
as it often provides interesting insights into fluid turbulence. In ID the Burgers equation is 

d t v + vd x v = vX7 2 v + f, (12) 

where / is the external force and the velocity v can have shocks since the system is 
compressible. In the unforced, inviscid case the Burgers equation has infinitely many 
conserved quantities, namely, J v n dx for all integers n. In the limit v — > we can 
use the Cole-Hopf transformation, v = d x ^, / = —d x F, and <3> = 2v In 6, to obtain 
d t Q = vd 2 Q + FQ/(2v), a linear partial differential equation (PDE) that can be solved 
explicitly in the absence of any boundary [38,39]. 

Passive scalars such as pollutants can be advected by fluids. These flows are governed 
by the advection-diffusion equation 

d t e + u.ve = K v 2 e + fg, (13) 

where 9 is the passive-scalar field, the advecting velocity field u satisfies the NS equation 7, 
and ig is an external force. The field 9 is passive because it does not act on or modify u. 
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Note that Eq.( 13) is linear in 6. It is possible, therefore, to make considerable analyti- 
cal progress in understanding the statistical properties of passive-scalar turbulence for the 
simplified model of passive-scalar advection due to Kraichnan [34,97]; in this model each 
component of is a zero-mean Gaussian random variable that is white in time; further- 
more, each component of u is taken to be a zero-mean Gaussian random variable that is 
white in time and which has the covariance 

(« i (x,t)« J -(x + r,t / )) = 2Dij6(t - t'); (14) 

the Fourier transform of Dij has the form 

^oc^ + ^-^V^^-^]; (15) 

q is the wave vector, L is the characteristic large length scale, r\ is the dissipation scale, 
and £ is a parameter. In the limit of L — > 00 and 77 — > we have, in real space, 

D ij (r)=D°5 ij -±d ij (r) (16) 

with 

d ij =D 1 r£[{d-l + Q5 iJ -^]. (17) 

D\ is a normalization constant and £ a parameter; for < £ < 2 equal-time passive-scalar 
structure functions show multiscaling [34]. 

We turn now to an example of a model for non-Newtonian flows. This model combines 
the NS equation for a fluid with the finitely extensible nonlinear elastic Peterlin (FENE-P) 
model for polymers; it is used inter alia to study the effects of polymer additives on fluid 
turbulence. This model is defined by the following equations: 

8 t u + (u.V)u = i/V 2 u + A V .[/(r P )C] - Vp; (18) 

TP 

d t C + u.VC = C.(Vu) + (Vu) T .C - f( rp ^ C - 1 . (19) 

Tp 

Here v is the kinematic viscosity of the fluid, /i the viscosity parameter for the solute 
(FENE-P), Tp the polymer relaxation time, p the solvent density, p the pressure, (Vu) T 
the transpose of (Vu), C a p = {R a Rp) the elements of the polymer-conformation tensor 
C (angular brackets indicate an average over polymer configurations), X the identity ten- 
sor with elements 6 a p, f(rp) = (L 2 - 3)/(L 2 - rf>) the FENE-P potential that ensures 
finite extensibility, rp = ^Tr(C) and L the length and the maximum possible extension, 
respectively, of the polymers, and c = + /x) a dimensionless measure of the polymer 
concentration [98]. 

The hydrodynamical partial differential equations (PDEs) discussed above are difficult 
to solve, even on computers via direct numerical simulation (DNS), if we want to resolve 
the large ranges of spatial and temporal scales that become relevant in turbulent flows. It 
is useful, therefore, to consider simplified models of turbulence that are numerically more 
tractable than these PDEs. Shell models are important examples of such simplified models; 
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they have proved to be useful testing grounds for the multiscaling properties of structure 
functions in turbulence. We will consider, as illustrative examples, the Gledzer-Ohkitani- 
Yamada (GOY) shell model [99] for fluid turbulence in three dimensions and a shell model 
for the advection-diffusion equation [100]. 

Shell models cannot be derived from the NS equation in any systematic way. They 
are formulated in a discretised Fourier space with logarithmically spaced wave vectors 
k n = fcoA™, A > 1, associated with shells n and dynamical variables that are the complex, 
scalar velocities u n . Note that k n is chosen to be a scalar: spherical symmetry is implicit 
in GOY-type shell models since their aim is to study homogeneous, isotropic turbulence. 
Given that k n and u n are scalars, shell models cannot describe vortical structures or enforce 
the incompressibility constraint. 

The temporal evolution of such a shell model is governed by a set of ordinary differen- 
tial equations that have the following features in common with the Fourier-space version 
of the NS equation [12]: they have a viscous-dissipation term of the form —vk\u n , they 
conserve the shell-model analogues of the energy and the helicity in the absence of viscos- 
ity and forcing, and they have nonlinear terms of the form ik n u n u n i that couple velocities 
in different shells. In the NS equation all Fourier modes of the velocity affect each other 
directly but in most shell models nonlinear terms limit direct interactions to shell velocities 
in nearest- and next-nearest-neighbour shells; thus direct sweeping effects, i.e., the advec- 
tion of the largest eddies by the the smallest eddies, are present in the NS equation but not 
in most shell models. This is why the latter are occasionally viewed as a highly simplified, 
quasi-Lagrangian representation (see below) of the NS equation. 

The GOY-model evolution equations have the form 

+ vk^u-a = i(a n u n+ iu n+2 b n u n -iu n+ i + c„u„_iu„_ 2 )* + /„, (20) 

where complex conjugation is denoted by *, the coefficients are chosen to be a n = k n , 
b n = —5k n -i, c n = — (1 — S)k n - 2 to conserve the shell-model analogues of the energy 
and the helicity in the inviscid, unforced case; in any practical calculation 1 < n < N, 
where N is the total number of shells and we use the boundary conditions u n = V n < 1 
or Vn > N; as mentioned above k n = A"fc and many groups use A = 2, S = 1/2, 
fco = 1/16, and N = 22. The logarithmic discretisation here allows us to reach very high 
Reynolds number, in numerical simulations of this model, even with such a moderate value 
of N. For studies of decaying turbulence we set /„ = 0, Vn; in the case of statistically 
steady, forced turbulence [45] it is convenient to use /„ = (1 + i)5 x 10~ 3 . For such a 
shell model the analogue of a velocity structure function is S p (k n ) — (\u(k n )\ p ) and the 
energy spectrum is E(k n ) = \u{k n )\ 2 /k n . 

It is possible to construct other shell models, by using arguments similar to the ones we 
have just discussed, for other PDEs such as the advection-diffusion equation. Its shell- 
model version is 




{9 n+2 u n+ i + n+ iu n+2 )] 



6 n -iu n +i) — 



(21) 
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For this model, the advecting velocity field can either be obtained from the numerical 
solution of a fluid shell model, like the GOY model above, or by using a shell-model ver- 
sion of the type of stochastic velocity field introduced in the Kraichnan model for passive- 
scalar advection [46]. A shell-model analogue for the FENE-P model of fluid turbulence 
with polymer additives may be found in Ref. [101]. 



3.1 Eulerian, Lagrangian, Quasi-Lagrangian frameworks 

The Navier-Stokes Eq.( 7) is written in terms of the Eulerian velocity u at position x and 
time t; i.e., in the Eulerian case we use a frame of reference that is fixed with respect to the 
fluid; this frame can be used for any flow property or field. The Lagrangian framework [5] 
uses a complementary point of view in which we fix a frame of reference to a fluid particle; 
this fictitious particle moves with the flow and its path is known as a Lagrangian trajectory. 
Each Lagrangian particle is characterised by its position vector ro at time to; its trajectory 
at some later time t is R = R(t; r , to) and the associated Lagrangian velocity is 



fIR. 

dt 



(22) 



We will also employ the quasi-Lagrangian [102,103] framework that uses the following 
transformation for an Eulerian field ip(r, t): 

^(r,i)=#- + R(i;r o ,0),t]; (23) 

here ip is the quasi-Lagrangian field and R(i; r , 0) is the position at time t of a Lagrangian 
particle that was at point r at time t = 0. 

As we have mentioned above, sweeping effects are present when we use Eulerian veloc- 
ities. However, since Lagrangian particles move with the flow, such effects are not present 
in Lagrangian and quasi-Lagrangian frameworks. In experiments neutrally buoyant tracer 
particles are used to obtain Lagrangian trajectories that can be used to obtain statistical 
properties of Lagrangian particles. 



4. Homogeneous Isotropic Turbulence: Phenomenology 

In 1941 Kolmogorov [56] developed his classic phenomenological approach to turbulence 
that is often referred to as K41. He used the idea of the Richardson cascade to provide an 
intuitive, though not rigorous, understanding of the power-law behaviours we have men- 
tioned in Sec. 2. We give a brief introduction to K41 phenomenology and related ideas; 
for a detailed discussion the reader should consult Ref. [33]. 

First we must recognise that there are two important length scales: (a) The large integral 
length scale L that is comparable to the system size and at which energy injection takes 
place; flow at this scale depends on the details of the system and the way in which energy 
is injected into it; (b) and the small dissipation length scale -qa below which energy dissi- 
pation becomes significant. The inertial range of scales, in which structure functions and 
energy spectra assume the power-law behaviours mentioned above (Sec. 2), lie in between 
L and rj; as Re increases so does the extent of the inertial range. 
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In K41 Kolmogorov made the following assumptions: (a) Fully developed 3D turbu- 
lence is homogeneous and isotropic at small length scales and far away from boundaries, 
(b) In the statistical steady state, the energy dissipation rate per unit volume e remains fi- 
nite and positive even when Re — > oo (the dissipative anomaly mentioned above), (c) A 
Richardson-type cascade is set up in which energy is transferred by the breakdown of the 
largest eddies, created by inherent instabilities of the flow, to smaller ones, which decay 
in turn into even smaller eddies, and so on till the sizes of the eddies become comparable 
to rjd where their energy can then be degraded by viscous dissipation. As Re — > oo all 
inertial-range statistical properties are uniquely and universally determined by the scale r 
and e and are independent of L, v and 77^. 

Dimensional analysis then yields the scaling form of the order-p structure function 

Sp A1 {r) « Ce p /V/ 3 , (24) 

since e has dimensions of (length) 2 (time)~ 3 . [It is implicit here that the eddies, at any 
given level of the Richardson cascade, are space filling; if not, e is intermittent and scale 
dependent as we discuss in Sec. 5 on multiscaling.] Thus (p 41 — p/3; for p = 2 
we get S^ 41 (r) ~ r 2 / 3 whose Fourier transform is related to the K41 energy spectrum 
E(k) K41 ~ fc~ 5 / 3 (left panel of Fig. 1). 

The prediction C;f 41 = 1, unlike all others K41 results, can be derived exactly for the 
NS equation in the limit Re — ► 00. In particular, it can be shown that [33,44] 

S 3 (£) » (25) 

an important result, since it is both exact and nontrivial. 

It is often useful to discuss K41 phenomenology by introducing ve, the velocity associ- 
ated with the inertial-range length scale £; clearly 

vt~t 1/3 l 1/3 . (26) 

The time scale tg ~ ^, the typical time required for the transfer of energy from scales of 
order I to smaller ones. This yields the rate of energy transfer 

n~£~f (27) 

Given the assumptions of K41, there is neither direct energy injection nor molecular dissi- 
pation in the inertial range. Therefore, the energy flux II becomes independent of £ and is 
equal to the mean energy dissipation rate e, which can now be written as 

e ~ vl/t (28) 

A similar prediction, for the two-point correlations of a passive-scalar advected by a 
turbulent fluid is due to Obukhov and Corsin; we shall not discuss it here but refer the 
reader to Ref. [104,105]. 

As we have mentioned above, the cascade of energy in 3D turbulence is replaced in 
2D turbulence by a dual cascade: an inverse cascade of energy from the injection scale 
to larger length scales and a forward cascade of enstrophy [35,36,74,75]. In the inverse 
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log k > p 

Figure 1. (Color online) (a) A representative log-log plot of the energy spectrum 
E(k) versus k, from a numerical simulation of the GOY shell model with 22 shells. 
The straight black line is a guide to the eye indicating K41 scaling fc~ 5//3 . (b) A plot of 
the equal-time scaling exponents £ p versus p, with error bars, obtained from the GOY 
shell model. The straight black line (color online) indicates K41 scaling p/3. 



cascade the energy accumulation at large length scales is controlled eventually by Ekman 
friction. The analogue of K41 phenomenology for this case is based upon physical argu- 
ments due to Rraichnan, Leith and Batchelor [75]. Given that there is energy injection 
at some intermediate length scale, kinetic energy get redistributed from such intermediate 
scales to the largest length scale. The scaling result for the two cascades gives us a kinetic 
energy spectrum that has a fc~ 5 / 3 form in the inverse-cascade inertial range and a fc~ 3 form 
(in the absence of Ekman friction) in the forward-cascade inertial range. 



5. From scaling to multiscaling 

In equilibrium statistical mechanics, equal-time and time-dependent correlation functions, 
in the vicinity of a critical point, display scaling properties that are well understood. For 
example, for a spin system in d dimensions close to its critical point, the scaling forms of 
the equal-time correlation function g(r; i, h) and its Fourier transform g(k; i, h), for a pair 
of spins separated by a distance r, are as follows: 



g(k;t,h)n V 1 fc2 y '-. (30) 

Here the reduced temperature t = (T — T c )/T c , where T and T c are, respectively, the 
temperature and the critical temperature, and the reduced field h — H/kgT,., with H 
the external field and ks the Boltzmann constant. The equal-time critical exponents fj, v 
and A are universal for a given universality class (the unconventional overbars are used 
to distinguish these exponents from the kinematic viscosity, etc.). The scaling functions 
G and G can be made universal too if two scale factors are taken into account [106]. 
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Precisely at the critical point (i = 0, h = 0) these scaling forms lead to power-law decays 
of correlation functions; and, as the critical point is approached, the correlation length £ 
diverges [e.g., as £ ~ t^~ v ^ if h = 0]. Time-dependent correlation functions also display 
scaling behaviour; e.g., the frequency (w) dependent correlation function has the scaling 
form to Eq. (30). 

g(k,w;t,h)& ^— . (31) 

This scaling behaviour is associated with the divergence of the relaxation time 

r~e, (32) 



referred to as critical slowing down; here z is the dynamic scaling exponent. 

In most critical phenomena in equilibrium statistical mechanics we obtain the simple 
scaling forms summarised in the previous paragraph. The inertial-range behaviours of 
structure functions in turbulence (Sees. 2 and 3) are similar to the power-law forms of 
these critical-point correlation functions. This similarity is especially strong at the level of 
K41 scaling (Sec. 4); however, as we have mentioned earlier, experimental and numerical 
work suggests significant multiscaling corrections to K41 scaling with the equal-time mul- 
tiscaling exponents ( p ^ ( p 41 ; in brief, multiscaling implies that ( p is not a linear function 
p; indeed [33] it is a monotone increasing nonlinear function of p (see right panel of Fig. 
1). The multiscaling of equal-time structure functions seems to be a common property of 
various forms of turbulence, e.g., 3D turbulence and passive-scalar turbulence. 

The multifractal model [33,107,108] provides a way of rationalising multiscaling cor- 
rections to K41. First we must give up the K41 assumption of only one relevant length 
scale I and the simple scaling form of Eq.( 28). Thus we write the equal-time structure 
function as 

S p (e) = C p (eer^(^-) s -, (33) 

where 8 P = C, P — p/3 is the anomalous part of the scaling exponent. We start with the 
assumption that the turbulent flow possesses a range of scaling exponents h in the set 
I = (h min ,h m ax)- For each h in this range, there is a set E/, (in real space) of fractal 
dimension D{h), such that, as £/L — > 0, 

8v e (r)~l h , (34) 

if r G Eft. The exponents (h min , h max ) are postulated to be independent of the mechanism 
responsible for the turbulence. Hence 

S p {i) ~ J dn(h)(e/L)P h+3 - D ^, (35) 

where the ph term comes from p factors of (l/L) in Eq. (34) and the 3 - D(h) term 
comes from an additional factor of (£/L) 3 ~ D ( h \ which is the probability of being within 
a distance of ~ £ of the set E^ of dimension D(h) that is embedded in three dimensions. 
The co-dimension D(h) and the exponents h min and h max are assumed to be universal 
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[33]. The measure dfi(h) gives the weight of the different exponents. In the limit l/L — > 
the method of steepest descent yields 

C P = inf h \ph + 3-D(h)]. (36) 

The K41 result follows from Eq. (36) if we allow for only one value of h, namely, h = 1/3 
and set D(h) = 3. For more details we refer the reader to [33,107,108]; the extension to 
time-dependent structure functions is given in Refs. [45,46,109]. 

Exact results for multiscaling can be obtained for the Kraichnan model of passive-scalar 
turbulence. We outline the essential steps below; details may be found in Ref. [34]. 

The second-order correlation function is defined as 



C 2 (l,t) = <0(x,t)0(x + l,t)>. (37) 

Here the angular brackets denote averaging over the statistics of the velocity and the force 
which are assumed to be independent of one another [34]. This equation of motion 

8 t C 2 {\, t) = (9 t 0(x, i)0(x + 1, t)> + (0(x, t)dtO(x + 1, t)> (38) 

is easy to solve by first by using the advection-diffusion equation and then using Gaussian 
averages to obtain [34] 

d t C 2 (l) = Dj-'dtKd - l)l d - 1+i C 2 (l)} + 2 K l 1 - d d l [l d - 1 d l C 2 (l)} + $(-!-), 

(39) 

where is the spatial correlation of the force [34] (notice that we now work with just 

the scalar I for the isotropic case). In the stationary state the time derivative vanishes on 
the left hand side. We impose the boundary conditions that, as I — > oo, C 2 (l) = 0, and 
C 2 (l) remains finite when I — > 0, whence 

i r°° r i-d r r r 

In the limit Id « I « L\, the second-order structure function has the following scaling 
form, 

S 2 (l) = 2[C 2 (0) - C 2 (l)} « {2 _Q( d _ 1)D mi 2 - £ , (41) 

i.e., equal-time exponents £f = 2 — ^; this result follows from dimensional arguments as 
well. For order-p correlation functions the equivalent of Eq. (38) can be written symboli- 
cally as [34] 

d t C p = -M P C P + D P C P + F® C p - 2 (42) 

where the operator M p is determined by the advection term, D p is the dissipative operator, 
and F is the spatial correlator of the force. In the limit of vanishing diffusivity, and in 
stationary state, the above equation reduces to 
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M P C P = F ® C p -2- (43) 

The associated homogeneous and inhomogeneous equations can be solved separately. By 
assuming scaling behaviour, we can extract the scaling exponent from simple dimensional 
analysis (superscript dim) to obtain 

Cf"-f(2-£). (44) 

The solution Z p (Xri, Ar 2 ...Ar p ) of the homogeneous part of Eq. (43) are called the 
zero-mode of the operator M p . The zero-modes have the scaling property 

Zp(Ari,Ar 2 ...Ar p ) - A c " r< %(ri, r 2 ...r p ). (45) 

Their scaling exponents Q ero cannot be determined from dimensional arguments. The 
exponents £p 6 ™ are a ls° called anomalous exponents. And for a particular order-p the 
actual scaling exponent is 

C P = min(C p ze ™, Cf") (46) 

This is how multiscaling arises in Kraichnan model of passive-scalar advection. The prin- 
cipal difficulty lies in solving the problem with a particular boundary condition. In recent 
times the following results have been obtained: Although the scaling exponents for the 
zero-modes has not been obtained exactly for any p, except for p = 2 (in which case 
the anomalous exponent is actually subdominant), perturbative methods have yielded the 
anomalous exponents. Also, it has been shown that the multiscaling disappears for £ > 2 
or £ < and that, although the scaling exponents are universal, the amplitudes depend on 
the force correlator and hence the structure functions themselves are not universal. These 
results have been well supported by numerical simulations. 

Several studies of the multiscaling of equal-time structure functions have been carried 
out as outlined above. By contrast there are fewer studies of the multiscaling of time- 
dependent structure functions. We give an illustrative example for the Kraichnan model 
of passive-scalar advection. For simplicity, we look at the Eulerian second-order time- 
dependent structure function which is defined, in Fourier space, as [46,1 10] 

j*(k,to,t) = <0(-k,t o )0(k,i)>. (47) 

In order to arrive at a scaling form for ^"(k, to,t), we look at its equation of motion: 

^ = (»H,<.)«). («) 

A spatial Fourier transform of the advection-diffusion equation (13) yields 

gg(k) 
dt 

so (48) maybe expanded as 

(tf*(k,t ,t) 



i J k jUj (q)e(k - q)d d q - /cfy/^k), (49) 



dt 



ikj J ((9(-k, t )uj(q)9(k - q, t))d d q - Kkjkj{9(-\si, h)9(k, t)). 



(50) 
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The above equation is solved with the help of Gaussian averaging. The first term reduces 
to 

f°° _ A 

(0(-]MoK(q)0(k-q,i)>= jf („ J .(i)^(t'))(0(-k,t o )^^0(k- q,^ 

(51) 

Equations (14) and (49) yield 
d^(k,t ,t) 



-2kikj / D l3 d d qF{^t^t). (52) 
Jo 



Since 2 f °° Dijd d q — D a (L) ~ L^, the equation of motion of the second-order structure 
function for the Eulerian field becomes 

dT°(r,t ,t) = L£ d 2 T (r,t o ,t) 
dt " dr 2 

whence [46] 

^(k,io,t) = ^,io)e- fc2i ' 4 . (54) 

Thus it is clear that within the Eulerian framework we get a simple dynamic scaling expo- 
nent 2 = 2. 

A similar analysis for the quasi-Lagrangian time-dependent structure function [46] gives 



dT(r,t ,t) _ ^^(r, t ,t) d.F(r, t ,t) 

^ ^ 9r . ar . dr . dr . 



A Fourier transform of Eq. (55) yields .F(k, to, i) oc exp[— t/r], wherer = fc^~ 2 , which 
implies a simple dynamic scaling exponent z = 2 — £ in the quasi-Lagrangian framework. 
In Sec. 6.2 we discuss dynamic scaling and multiscaling in shell models. 



6. Numerical Simulations 

Numerical studies of the models described in Sec. 3 have contributed greatly to our un- 
derstanding of turbulence. In this Section we give illustrative numerical studies of the 3D 
Navier-Stokes equation (Sec. 6.1), GOY and advection-diffusion shell models (Sec. 6.2), 
the 2D Navier-Stokes equation (Sec. 6.3), the ID Burgers equation (Sec. 6.4) and the 
FENE-P model for polymer additives in a fluid (Sec. 6.5). 



6.1 3D Navier-Stokes Turbulence 

We concentrate on the statistical properties of homogeneous, isotropic turbulence, so we 
restrict ourselves periodic boundary conditions. Even with these simple boundary condi- 
tions, simulating these flows is a challenging task as a wide range of length scales has to be 
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resolved. Therefore, state-of-the-art numerical simulations use pseudo-spectral methods 
that solve the Navier-Stokes equations via Fast Fourier transforms [111,112] typically on 
supercomputers. For a discussion on the implementation of the pseudo-spectral method 
we refer the reader to Refs. [111,112]. We outline this method below: (a) Time marching 
is done by using either a second-order, slaved Adams-Bashforth or a Runge-Kutta scheme 
[113]. (b) In Fourier space the contribution of the viscous term is -vk 2 u. (c) To avoid 
the computational costs of evaluating the convolution because of the non-linear term, it is 
first calculated in real space and then Fourier transformed; hence the name pseudo-spectral 
method, (d) In Fourier space the discretized Navier-Stokes time evolution us 

u n+1 = cxp(-vk 2 5t)u n + 1 ~ exp(-z/fc 2 Jt) _ (1/2)A p-i] 

vkt 

where n is the iteration number, Af indicates the non-linear term, and = (Sij—kikj/k 2 ) 
is the transverse projector which guarantees incompressibility. (e) To suppress aliasing 
errors we use a 2/3 dealiasing scheme [112]. 

We give illustrative results from a direct numerical simulation DNS with 1024 3 that we 
have carried out. This study uses the stochastic forcing of [1 14] and has attained a Taylor 
microscale Reynolds number Re\ ~ 100, where Re\ — u rms \/v\ u rms = ^/2E/3 
is the root-mean-square velocity and the Taylor microscale A = E(k) I k 2 E(k). 
For state-of-the-art simulations with up to 4096 3 collocation points we refer the reader to 
Ref. [79]. As we had mentioned in Sec. 2, regions of high vorticity are organised into 
slender tubes. These can be visualised by looking at isosurfaces of \ui\ as shown in the 
representative plots of Figs 2 and 3. The right panel of Fig. 2 shows the PDF of \u>\; 
this has a distinctly non-Gaussian tail. The structure of high-|cj| vorticity tubes shows 
up especially clearly in the plots of Fig. 3, the second and third panels of which show 
successively magnified images of the central part of the first panel (for a 4096 3 version see 
Ref. [79]). 




One method to look at these structures is to study the joint PDF of the invariants Q = 
—tr(A 2 )/2 and R = — tr(A 3 )/3 of the velocity gradient tensor. The zero-discriminant or 
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Figure 3. (Color online) (Left) Isosurface plot of |ojj with \u\ equal to one standard 
deviation more than its mean value. (Center) A magnified version of the central part of 
the panel on the left. (Right) A magnified version of the central part of the panel in the 
middle. 
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Figure 4. (Color online) (Left) Joint PDF P(Q*,R*) of 7?* = H/is^s^) 3 ^ and 
Q* — Q/( s ij s ij) calculated from our DNS. The black curve represents the zero-dis- 
criminant (or Vieillefosse) line 27i? 2 /4 + Q [i — 0. The contour levels are logarith- 
mically spaced. (Right) PDF of the a>component of the velocity (here a denotes the 
standard deviation); the parabolic curve is a Gaussian that is drawn for comparison. 
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Figure 5. (Color online) (Left) The compensated energy spectrum k 5 ^ 3 E(k) versus 
krj, where r\ is the dissipation scale from our DNS (see text). (Right) PDFs of velocity 
increments that show marked deviation from Gaussian behaviour (innermost curve), 
especially at small length scales; the outermost PDF is for the velocity increment with 
the shorter length scale. 



Vieillefosse line D = 27i? 2 /4 + Q 3 = divides the QR plane in different regions. The 
region with D > is vorticity dominant (one of the eigenvalues of A is greater than zero 
whereas the other two eigenvalues are imaginary); the region D < is strain dominated 
(all the eigenvalues of A are real). The regions D > and D < can be further divided 
into two more quadrants depending upon the sign of the eigenvalues. In the left panel of 
Fig. 4 we show a representative contour plot of the joint PDF P(Q* , R*) obtained from 
our DNS. The shape of the contour is like a tear-drop, as in experiments [63], with a tail 
along the line D = in the region where R* > and Q* < 0. The plot indicates that, in a 
numerical simulation, most of the structures are vortical but there also exist regions of large 
strain. For a more detailed discussion of the above classification of different structures we 
refer the reader to [63,1 15]. 

The left panel of Fig. 5 shows a plot of the compensated energy spectrum k 5 ^ 3 E(k) 
versus krj (77 is the dissipation scale in our DNS). The flat portion at low kr\ indicates 
agreement with the K41 form E K41 (k) ~ fc -5 / 3 . There is a slight bump after that; this 
is referred to as a bottleneck (see Ref. [116] and Sec 6.4); the spectrum then falls in the 
dissipation range. The right panel of Fig. 5 shows PDFs of velocity increments at different 
scales r. The innermost curve is a Gaussian for comparison; the non-Gaussian deviations 
increase as r decreases. 

We do not provide data for the multiscaling of velocity structure functions in the 3D 
Navier-Stokes equation. We refer the reader to Ref. [60] for a recent discussion of such 
multiscaling. Often the inertial range is quite limited in such studies. This range can be 
extended somewhat by using the extended-self-similarity (ESS) procedure [117] in which 
the slope of a log-log plots of the structure function S p versus S q yields the exponent ratio 
Cp/(q', this procedure is especially useful if q — 3 since (3 = 1 for the 3D Navier-Stokes 
case. We illustrate the use of this ESS procedure in Sec. (6.3) on 2D turbulence. 

The methods of statistical field theory have been used with some success to study the 
statistical properties of a randomly forced Navier-Stokes equation [25,26,30,31]. The 
stochastic force here acts at all length scales; it is Gaussian and has a Fourier-space co- 
variance proportional to fc 1_y . For y > 0, a simple perturbation theory leads to infrared 
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divergences; these can be controlled by a dynamical renormalization group for sufficiently 
small y; for y = 4 this yields a K41-type fc~ 5 / 3 spectrum at the one-loop level. This value 
of y is too large to trust a low-y, one-loop result; also, for y > 3, the sweeping effect 
leads to another singularity [118]. Nevertheless, this randomly forced model has played an 
important role historically. Thus it has been studied numerically via the pseudo-spectral 
method [119,120]. These studies have shown that, even though the stochastic forcing de- 
stroys the vorticity tubes that we have described above, it yields multiscaling of velocity 
structure that is consistent, for y = 4, with the analogous multiscaling in the conventional 
3D Navier-Stokes equation, barring logarithmic corrections. We will discuss the analogue 
of this problem for the stochastically forced Burgers equation in Sec. 6.4. 



6.2 Shell Models 

Even though shell models are far simpler than their parent partial differential equations 
(PDEs), they cannot be solved analytically. The multiscaling of equal-time structure func- 
tions in such models has been investigated numerically by several groups; an overview of 
earlier work and details about numerical methods for the stiff shell-model equations can 
be found in Refs. [45,46,121]. An illustrative plot of equal-time multiscaling exponents 
for the GOY shell model is given in the right panel of Fig. 1 . 

We devote the rest of this Subsection to a discussion of the dynamic multiscaling of 
time-dependent shell-model structure functions that has been elucidated recently by our 
group [45,46,109,1 10]. So far, detailed numerical studies of such dynamic multiscaling has 
been possible only in shell models. We concentrate on time-dependent velocity structure 
functions in the GOY model and their passive-scalar analogues in the advection-diffusion 
shell model. 

In a typical decaying-turbulence experiment or simulation, energy is injected into the 
system at large length scales (small fc); it then cascades to small length scales (large k); 
eventually viscous losses set in when the energy reaches the dissipation scale. We will 
refer to this as cascade completion. Energy spectra and structure functions show power- 
law forms like their counterparts in statistically steady turbulence. It turns out [46] that 
the multiscaling exponents for both equal-time and time-dependent structure functions are 
universal in so far as they are independent of whether they are measured in decaying tur- 
bulence or the forced case in which we get statistically steady turbulence. 

Furthermore, the distinction between Eulerian and Lagrangian frameworks assumes spe- 
cial importance in the study of dynamic multiscaling of time-dependent structure functions. 
Eulerian-velocity structure functions are dominated by the sweeping effect that lies at the 
heart of Taylor's frozen-flow hypothesis; this relates spatial and temporal separations lin- 
early (see Sec. 2) whence we obtain trivial dynamic scaling with dynamic exponents 
Zp = 1 for all p, where the superscript £ stands for Eulerian. By contrast, we expect 
nontrivial dynamic multiscaling in Lagrangian or quasi-Lagrangian measurements. Such 
measurements are daunting in both experiments and direct numerical simulations; how- 
ever, they are possible in shell models. As we have mentioned in Sec. 3, shell models 
have a quasi-Lagrangian character since they do not have direct sweeping effects. Thus we 
expect nontrivial dynamic multiscaling of time-dependent structure functions in them. 

Indeed, we find that [45,46,103] that, given a time-dependent structure function, we can 
extract an infinity of time scales from it. Dynamic scaling Ansatze [cf., Eq. (4)] can then 
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be used to extract dynamic multiscaling exponents. A generalisation of the multifractal 
model then suggests linear relations, referred to as bridge relations, between these dynamic 
multiscaling exponents and their equal-time counterparts. These can be related to equal- 
time exponents via bridge relations. We show how to check these bridge relations in shell 
models. However, before we present details, we must define time-dependent structure 
functions precisely. 

The order-p, time-dependent, structure functions, for longitudinal velocity increments, 
Suu (x, r, t) = [u(x + r, t) — u(x, t)] and passive-scalar increments, <5#(x, t, r) = #(x + 
r, t) — 6*(x, t) are defined as 



T p {r, t p }) = ([5u\i (x, h,r)... Su\\ (x, t p , r)]) 



(56) 



and 



tf(r,t u ...,t p ) =< [50(x.,t 1 ,r)...86{x,t p ,T)] >; 



(57) 



i.e., fluctuations are probed over a length scale r which lies in the inertial range. For 
simplicity, we consider t\ — t and t% = . ■ . = t p = in both Eq. (56) and Eq. (57). 
Given T u (r,t) and J- 6 (r, t), we can define the order-p, degree-M, integral-time scales 
and derivative-time scales as follows [46]: 



1 r°° 



(l/M) 



(58) 



W Jo 



T e p (r,t)t^ M ^dt 



(l/M) 



(59) 



1 d M T^{r,t) 
S%{r) dt M 



(-l/M) 



(60) 



P,M 



(r,t) 



1 d M ^(r,t) 



S*(r) dt M 



(-1/M) 



(61) 



Integral-time dynamic multiscaling exponents z pM for fluid turbulence can be defined 



r z P ,M an( j me derivative-time ones z^ M by T^jfi (r,t) ~ r z p> M . They 



satisfy the following bridge relations [46]: 

1 + [Cp-m - Q/M; 



*P,M 



I,u 

pM 



(62) 



pM 



1 + [C P - ( p +m}/M. 



(63) 



For passive-scalars advected by a turbulent velocity field, the corresponding dynamic mul- 

10 1,0 DO D '° 

tiscaling exponents are defined as T pM (r, t) oc r z p> M and T p M (r, t) oc r z p- M ; they satisfy 
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the following bridge relations involving the scaling exponents £m of equal-time, order- M 
structure functions of the advecting velocity field: 



J ,(7 -i 

Z p,M — 1 



Cm 
M' 



D.e _ , _ C-m 

Z VM - 1 M ■ 



(64) 



These bridge relations, unlike Eq. (62) and Eq. (63), are independent of p. [Recall that, 
for the Rraichnan model, we have already shown in Sec. 5 that we get simple dynamic 
scaling.] 

GOY-model equal-time structure functions and their associated inertial-range exponents 
are defined as follows: 



^(*») = ([«n(*K(*)] P/a )~*n C '- 
The time-dependent structure function are 

F?[k n ,t , t) = ([u n (t )u*(t + i)f /2 



(65) 



(66) 



We evaluate these numerically for the GOY shell model [numerical details may be found 
in Refs. [45,46]], extract integral and derivative time scales from them and thence the 



exponents z^ and z^^", respectively, from slopes of log-log plots of Tp'"(n) versus k n 



I,u, 



(right panel of Fig. 6) and of X^" (n) versus k n (right panel of Fig. 7). 




t log k 

Figure 6. (Color online) (a) A representative plot of the normalised fourth order 
time-dependent structure function versus the dimensionless time r obtained from the 
GOY shell model. The plots are for shells 4, 6, and 8 (from top to bottom), (b) A log-log 
plot of T\ '"(n) versus k (for convenience, we have dropped the subscript n in the label 
of the x-axis in the figure); a linear fit gives the dynamic mulstiscaling exponent 2 4 '". 

There is excellent agreement (within error bars) of the multiscaling exponents z^\ and 



,", obtained from our simulations, with the values computed from the appropriate bridge 



D. 
Z p,2 

relations using the equal-time exponents, ( p . 

For the passive-scalar case, the equal-time order-p structure functions is 



s e p (k n ) = ([6(t)e* n (t)} p/2 



-c 



(67) 
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Figure 7. (Color online) (a) A representative plot of the normalised sixth order 
time-dependent structure function versus the dimensionless time r obtained from the 
GOY shell model. The plots are for shells 4, 6, and 8 (from top to bottom), (b) A 
log-log plot of T®2 U ( n ) versus k (for convenience, we have dropped the subscript n in 
the label of the x-axis in the figure); a linear fit gives the dynamic multiscaling exponent 

D,u 
Z 6,2 ■ 

and its time-dependent version is 

F*(kn, t Q , t) =< [e n (t )e*(t + t)\p/ 2 > . (68) 

We consider decaying turbulence here with to a time origin. It is useful now to work 
with the normalised time-dependent structure function, Qp(n, to,i) = pa^' t ° ) • F° r ^ s 
case of passive-scalars advected by a velocity field which is turbulent (a solution of the 
GOY model), we calculate the integral (for M = 1) and derivative time scales (for M = 2) 
corresponding to Eq.(58) and Eq.(60), respectively. The slope of a log-log plot of T p \{n) 

— z J,e 

vs fc„ yields the integral time scale exponent, z p \, since T p \ (n) cx k n p ' 1 ■ Likewise, from 

plots of the derivative time scales we extract the exponent z p ' 2 . For a detailed discussion 
on dynamic multiscaling in this model we refer the reader to Refs. [46,109]. 

6.3 2D Navier-Stokes Turbulence 

We now consider illustrative numerical calculations for the 2D NS equations (9)-( 11). We 
begin with periodic boundary conditions for which we can use a pseudo-spectral method 
similar to the one given in the previous Subsection for the 3D NS case. We study decay- 
ing turbulence first with the source function / (the z component of the curl of some force 
V x F) set to 0. We use 1024 2 collocation points and the standard 2/3 dealiasing pro- 
cedure; for time marching we use a second-order Runge-Kutta scheme [113]. Our initial 
condition |w(/c)| 2 = fc _3 exp(— k 2 ) leads to a forward cascade. We seed the flow with 
Lagrangian tracers and use a cubic spline interpolation method to calculate their trajecto- 
ries [122]. Representative plots from our from our DNS are shown in Fig. 8. The first 
part (Fig. 8a) shows a compensated energy spectrum k 3 E(k) for the case with no Ekman 
friction. Figure 8b, from a DNS with Ekman friction cue = 0.1, Kolmogorov forcing [89], 
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and periodic boundary conditions, shows a trajectory of a Lagrangian tracer superimposed 
on a pseudocolour plot of the vorticity field at time t = 100; the tracer starts at the point 
marked with a circle (t = 0) and ends at the star (t = 100). For a state-of-the-art simulation 
that resolves both forward and inverse cascades in a forced DNS of 2D turbulence we refer 
the reader to Ref. [123]; such DNS studies have also investigated the scaling properties 
of structure functions and have provided some evidence for conformal invariance in the 
inverse cascade inertial range [124]. 




log k x 
Figure 8. (Color online) (a) A log-log plot of the compensated energy spectrum 
k i E{k) versus k from our DNS, of resolution 1024 2 , of two dimensional decaying 
turbulence with periodic boundary conditions. The flat region indicates a scaling form 
E(k) ~ k~ s . (b) The trajectory of a single Lagrangian particle over a time of order 100 
in a two-dimensional flow with drag and forcing. The starting point of the trajectory is 
in the middle of the box and is indicated by a red circle; the end point is indicated by a 
blue star. The trajectory is superimposed on a pseudocolor plot of the vorticity field cor- 
responding to the time at the end of the Lagrangian trajectory. The figure corresponds to 
a forced DNS of resolution 1024 2 with periodic boundary conditions, statistical steady 
state, and with a coefficient of Ekman friction ole = 0.1. 

We end with an illustrative example of a recent DNS study [89] that sheds light on 
the effect of the Ekman friction on the statistics of the forward cascade in wall-bounded 
flows that are directly relevant to laboratory soap-film experiments [125-128]. The de- 
tails of this DNS are given in Ref. [89]. In brief, uj is driven to a statistical steady state 
by a deterministic Kolmogorov forcing F u = ki n jFocos(ki n jX), with Fq the amplitude 
and ki n j the wavenumber on which the force acts; no-slip and no-penetration boundary 
conditions are imposed on the walls. The important non-dimensional control parame- 
ters are the Grashof number Q = 2ir\\F UJ \\2 / (kf n j pv 2 ) and the non-dimensional Ekman 
friction 7 = cce / (kf n jV), where we non-dimensionalize F u by 27r/(fcj nj -||.F a ,||2), with 
H-FLH2 = il a l-^wpdx) 1 / 2 and the length and time scales are made non-dimensional by 
scaling x by and t by k^/v. We use a fourth-order Runge-Kutta scheme for time 
marching and evaluate spatial derivatives via second-order and fourth-order, centered, fi- 
nite differences, respectively, for points adjacent to the walls and for points inside the 
domain. The Poisson equation is solved by using a fast-Poisson solver [113] and u is 
calculated at the boundaries by using Thorn's formula [89]. 

Since Kolmogorov forcing is inhomogeneous, we use the decomposition ip = (ip) + ip' 
and ui = (u>) + uo' , where the angular brackets denote a time average and the prime the 
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fluctuating part to calculate the order-p velocity and vorticity structure functions. Since 
this is a wall-bounded flow, it is important to extract the isotropic parts of these structure 
functions [89,129]. Furthermore, given our resolution (2049 2 ), it becomes necessary to use 
the ESS procedure to extract exponent ratios. Illustrative log-log ESS plots for velocity, 
S P (R), and vorticity, S£(R), structure functions are shown in the left and right panels, 
respectively, of Fig. 9; their slopes yield the exponent ratios that are plotted versus the order 
p in Fig. 10. The Rraichnan-Leith-Batchelor (KLB) predictions [75] for these exponent 
ratios, namely, C^ LS /Cf LB ~ r p/2 and Cp ' KLB /C% ' KLB ~ r°, agree with our values for 
Cp/C2 but not C,^ /(,%'■ velocity structure functions do not display multiscaling [left panel 
of Fig. 10] whereas their vorticity analogs do [note the curvature of the plot in the right 
panel of Fig. 10]. Similar results have been seen in DNS studies with periodic boundary 
conditions [130,123]. Additional results for PDFs of several properties can be obtained 
from our DNS [89]; these are in striking agreement with experimental results [126]. 




Figure 9. (Color online) (Left) Log-log ESS plots of the isotropic parts of the order-p 
velocity structure functions S P (R) versus fife (it); p = 3 (purple line with dots), p = 4 
(red line with square), p = 5 (green line with triangles), and p = 6 (blue line with 
circles). According to the KLB prediction S P (R) ~ Rp. (Right) Log-log ESS plots 
of the isotropic parts of the order-p vorticity structure functions S P (R) versus Sa(it); 
p = 3 (purple line with stars), p — 4 (red line with square), p = 5 (green line with 
triangles), and p — 6 (blue line with circles). 




1 2 3 4 5 6 



Figure 10. (Color online) (Left) Plots of the exponent ratios Q p / £2 versus p for the ve- 
locity differences. (Right) Plots of the exponent ratios / versus p for the vorticity 
differences. 
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6.4 The One dimensional Burgers Equation 



In this Subsection we present a few representative numerical studies of the ID Burgers 
equation. The first of these uses a pseudo-spectral method with 2 14 collocation points, 
the 2/3 dealising rule, and a fourth-order Runge-Kutta time-marching scheme. In the sec- 
ond study of a stochastically forced Burgers equation (see below) we use a fast-Legendre 
method that yields results in the zero-viscosity limit [131]. 

For the Burgers equation with no external forcing and sufficiently well-behaved initial 
conditions, the velocity field develops shocks, or jump discontinuities, which merge into 
each other with time. The time at which the first shock appears is usually denoted by t*. 
For all times greater than t*, it is possible to calculate, analytically, the scaling exponents 
£ p for the equal-time structure functions via S p = ([u(x + r, t) — u(x)] p ) ~ C p |r| p + Cp|r|, 
where the first term comes from the ramp, and the second term comes from the probability 
of having a shock in the interval \r\. As a consequence of this we have bifractal scaling : 
for < p < 1 the first term dominates leading to ( p = p and for p > 1 the second one 
dominates giving £ p = 1- This leads to an energy spectrum E(k) ~ fc~ 2 . Representative 
plots from our pseudo-spectral DNS, with v = 10 -3 and an initial condition u(x) =sin(a;) 
(for which t* = 1) are shown in Fig. 11; the left panel shows plots of the velocity field at 
times t = 0, 1, and £ = 1.5 and the right panel the energy spectrum at t = 1. 




log k 



Figure 11. (Color online) (Left) Snapshots of the solution of the Burgers equation 
obtained from our DNS with initial condition u(x) =sin x at times t = (blue), t = 1 
(black) and t = 2 (red). (Right) A representative log-log plot of E(k) versus k, at time 
t — 1 for the Burgers equation with initial conditions u{x) = sin x. 

The stochastically forced Burgers equation has played an important role in 
renormalization-group studies [131]. In particular, consider a Gaussian random force 
f(x, t) with zero mean and the following covariance in Fourier space: 

(f(k 1 ,t 1 )f{k 2 ,t2))=2D o \k\ 5(t 1 -t2)5(k 1 +k 2 y, (69) 

here f(k,t) is the spatial Fourier transform of f(x, t), Do is a constant, and the scaling 
properties of the forcing is governed by the exponent (3. For positive values of (3, the 
Burgers equation can be studied by using renormalization-group techniques; specifically, 
for (3 = 2 one recovers simple (Kardar-Parisi-Zhang or KPZ) scaling with the equal-time 
exponent ( p = p. It was hoped that forcing with negative values of (3 (in particular (3 = 
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— 1), which cannot be studied by renormalization-group methods, might yield multiscaling 
of velocity structure functions. 

However, our high-resolution study [131], which uses a fast-Legendre method, has 
shown that the apparent multiscaling of structure functions in this stochastic model might 
arise because of numerical artifacts. The general consensus is that this stochastically forced 
Burgers model should show bifractal scaling. In Fig. 12 we present representative plots 
of the velocity field (left panel, blue curve) and the scaling exponents (right panel) for this 
model. We have obtained the data for these figures by using a fast-Legendre method with 
2 18 collocation points. 




_ p 
Figure 12. (Color online) (Left) A snapshot of the velocity field (jagged line in blue) 
in steady state and the force in red from our fast-Legendre method DNS of the stochas- 
tically forced Burgers equation. (Right) A representative plot of the exponents ( p , with 
error-bars, for the equal-time velocity structure functions of the stochastically forced 
Burgers equation; bifractal scaling is shown by the black solid line; the deviations from 
this are believed to arise from artefacts (see text). 

Numerical studies of the Burgers equation have also proved useful in elucidating bot- 
tleneck structures in energy spectra [132,133](cf., the spectrum in the left panel of Fig. 
5). It turns out that such a bottleneck does not occur in the conventional Burgers equation. 
However, it does [134] occur in the hyperviscous one, in which usual Laplacian dissipa- 
tion operator is replaced by its a* power; this is known as hyperviscosity for a > 1. We 
show a representative compensated energy spectrum for the case a = 4 in the left panel of 
Fig. 13. We have obtained this from a pseudo-spectral DNS with 2 12 collocation points. 
The oi — ► oo limit is very interesting too since, in this limit, the hyperviscous Burgers 
equation maps on to the Galerkin-truncated version of the inviscid Burgers equation. In 
this Galerkin-truncated inviscid case, the Fourier modes thermalise [135,136]; in a com- 
pensated energy spectrum this shows up as E(k) ~ fc 2 , for large k [see the right panel of 
Fig. 13 for the case a — 200]. Such thermalisation effects in the Galerkin-truncated Euler 
equation have also attracted a lot of attention [137]; and the link between bottlenecks and 
thermalisation has been explored in our recent work [134] to which we refer the interested 
reader. 
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Figure 13. (Color online) (Left) A representative log-log plot of a bottleneck in the 
compensated energy spectrum k 2 E(k) of a hyperviscous Burgers equation with a — 4. 
(Right) A representative log-log plot of k 2 E(k) versus k for a — 200 at time t = 30. 
We see clear signatures of thermalization at large k (see text). 



6.5 Turbulence with Polymer Additives 

In this Subsection, we present a few results from our numerical study [138] of the analogue 
drag reduction by polymer additives in homogeneous, isotropic turbulence. This requires 
a DNS of considerably greater complexity than the ones we have described above. A 
naive pseudospectral method cannot be used for the FENE-P model given in Eqs. (18) 
and (19): the polymer conformation tensor C is symmetric and positive definite; however, 
in a practical implementation of the pseudo-spectral method it loses this property. We 
have employed a numerical technique that uses a Cholesky decomposition to overcome 
this problem; we refer the reader to Ref. [138] for these details. 

Our recent DNS of this model has shown that the natural analogue drag reduction in 
decaying, homogeneous, isotropic turbulence is dissipation reduction; the percentage re- 
duction DR can be defined as 



here the superscripts / and p stand, respectively, for the fluid without and with polymers 
and the superscript m indicates the time t m at which the cascade is completed. The de- 
pendence of DR on the polymer concentration parameter c and the Weissenberg number 
may be found in Ref. [138]. Here we show how the addition of polymers reduces small- 
scale structures in the turbulent flow: By a comparison of the isosurfaces of \ui\ in the 
left (without polymers) and right (with polymers) panels of Fig. 14, we see that slender 
vorticity filaments are suppressed by the polymers; this is in qualitative agreement with ex- 
periments [93]. The PDFs of | u>\, with and without polymers (left panel of Fig. 15) confirm 
that regions of large vorticity are reduced by polymers. The right panel of Fig. 15 shows 
how the polymers modify the energy spectrum in the dissipation range; this behaviour has 
been seen in recent experiments [96], which study the second-order structure function that 
is related simply to the energy spectrum. For a full discussion of these and related results 
we refer the re ader to Ref . [101,138]. 




(70) 
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Figure 14. (Color online) Constant- \cj\ isosurfaces for \ui\ = (|o>|) + a at cascade 
completion without and (Right) with polymers (c = 0.4); is the mean and a the 
standard deviation of |w|. 




Figure 15. (Color online) (Left) PDF of uo at cascade completion without (c = 0) 
and with polymers (c = 0.4). Note that regions of large vorticity are reduced on the 
addition of polymers. (Right) Representative plots of the energy spectra E p ' m {k) or 
E^' m (k) versus k for c = 0.1 (blue dashed line) and c = 0.4 (solid line). 
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7. Conclusions 



Turbulence provides us with a variety of challenging problems. We have tried to give 
an overview of some of these, especially those that deal with the statistical properties of 
turbulence. The choice of topics has been influenced, of course, by the areas in which 
we have carried out research. For complementary, recent overviews we refer the reader to 
Refs. [1-3]; we hope the other reviews and books that we have cited to will provide the 
reader with further details. 

We would like to thank CSIR, DST, and UGC (India) for support, and SERC (HSc) for 
computational resources. 
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